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Abstract 

Continuum elastic models that account for membrane thickness variations are especially useful in the 
description of nanoscale deformations due to the presence of membrane proteins with hydrophobic mis- 
match. We show that terms involving the gradient and the Laplacian of the area per lipid are significant 
and must be retained in the effective Hamiltonian of the membrane. We reanalyze recent numerical data, 
as well as experimental data on gramicidin channels, in light of our model. This analysis yields consistent 
results for the term stemming from the gradient of the area per molecule. The order of magnitude we find 
for the associated amplitude, namely 13-60 mN/m, is in good agreement with the 25 mN/m contribution 
of the interfacial tension between water and the hydrophobic part of the membrane. The presence of this 
term explains a systematic variation in previously published numerical data. 

Introduction 

As basic constituents of cell membranes, lipid bilayers [T] play an important role in biological processes, 
not as a passive background, but rather as a medium that responds to and influences, albeit in a subtle 
way, the behavior of other membrane components, such as membrane proteins [2]- The coupling between 
the lipid bilayer and guest molecules does not occur by the formation of chemical bonds, but rather by 
a deformation of the membrane in its entirety. To describe it, one must resort to concepts developed in 
soft matter physics for the understanding of self-assembled systems. 

At length scales much larger than their thickness, the elasticity of lipid bilayers is well described by 
the Helfrich model [3]. However, nanometer-sized inclusions, such as membrane proteins, deform the 
membrane over smaller length scales. In particular, some transmembrane proteins have a hydrophobic 
part with a thickness slightly different from that of the hydrophobic part of the membrane. Due to this 
hydrophobic mismatch, the hydrophobic core of the membrane locally deforms [1HS] • As this deformation 
affects the thickness of the membrane, and as its characteristic amplitude and decay length are both of 
a few nanometers [7], it cannot be described using the Helfrich model. In fact, since the range of such 
deformations is of the same order as membrane thickness, one can wonder to what extent continuum 
clastic models in general still apply, and what level of complexity is required for an accurate description. 
In particular, which terms must be retained in a deformation expansion of the effective Hamiltonian? 

Experimental data is available for the gramicidin channel [5] , a transmembrane protein formed by two 
protein monomers. The channel being large enough for the passage of monovalent cations, conductivity 
measurements [9] can detect its formation and lifetime, which are directly influenced by membrane prop- 
erties. The gramicidin channel can therefore act as a local probe for bilayer elasticity on sub-nanometer 
scales (see, e.g., Ref. [10]). Motivated by this opportunity, sustained theoretical investigations have been 
conducted in order to construct a model describing membrane thickness deformations [71H1ITT5] . Recently, 
detailed numerical simulations have been performed, giving access both to the material constants involved 
in elastic models and to the membrane shape close to a mismatched protein |14H16j . This numerical data 
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provides a good test for theoretical models. 

In this article, we put forward a modification to the models describing membrane thickness defor- 
mations. We argue that contributions involving the gradient (and the Laplacian) of the area per lipid 
should be accounted for in the effective Hamiltonian per lipid from which the effective Hamiltonian of 
the bilayer is constructed, following the approach of Refs. (THUS]- We show that these new terms cannot 
be neglected, as they contribute to important terms in the bilayer effective Hamiltonian. We discuss 
the differences between our model and the existing ones. We compare the predictions of our model with 
numerical data giving the profile of membrane thickness close to a mismatched protein [14H16] , and with 
experimental data on gramicidin lifetime |17j and formation rate |18j . 

Results: Membrane model 

We consider a bilayer membrane constituted of two identical monolayers, labeled by + and — , in contact 
with a reservoir of lipids with chemical potential [i. We write the effective Hamiltonian per molecule in 
monolayer ± as 

ft = \% ( S± - S o) 2 ±hH ± ± #(E± - + h {H±f 

+ f K K* + a (VS ± ) 2 + ft V 2 S ± + C (V 2 S ± ) 2 -n, (1) 

where is the area per lipid, while is the local mean curvature of the monolayer, and is its local 
Gaussian curvature (denoting by cf and cf the local principal curvatures [TDJ of the monolayer, we have 
i/ ± = (cf + cf)/2 and = cfcf). All these quantities are defined on the hydrophilic-hydrophobic 
interface of each monolayer. Eq. [T] corresponds to an expansion of for small deformations around 
the equilibrium state where the membrane is fiat and where each lipid has its equilibrium area £o- Any 
constant term in the free energy per lipid is included in a redefinition of the chemical potential fi. From 
now on, we will consider small deformations of an infinite flat membrane and we will work in the Monge 
gauge, so H ± ~ V 2 h ± /2 and K^ 1 ~ d 2 h ± dyh ± — (d x d y h ± ) 2 = det(d i djh ± ), where z = h ± (x, y) represents 
the height of the hydrophilic-hydrophobic interface of each monolayer with respect to a reference plane 
(x, y). The upper monolayer is labeled by + and the lower one by — . Many constants involved in Eq. [1] 
can be related to the constitutive constants of a monolayer: /q So = K a /2 is the compressibility modulus 
of the monolayer, /2/(2£q) = kq/2 is its bending rigidity, /jf/Eo = re/2 is its Gaussian bending rigidity, 
/1//2 = Co is its spontaneous (total) curvature, and /1//2 = c' Q is the modification of the spontaneous 
(total) curvature due to area variations (see Methods, Sec. II. ip . 

In the case where a = ft = ( = 0, Eq. |TJ is equivalent to the model of Rcf. [15], which is the basis 
of that developed in Refs. [I2HT5] . To our knowledge, existing membrane models including the area per 
lipid (or, equivalently, the two-dimensional lipid density) do not explicitly feature terms in the gradient, 
or Laplacian, of this variable [2D]. The possibility of an independent term proportional to the squared 
thickness gradient was however suggested on symmetry grounds in Ref. [2T] , while pointing that it could 
arise from the specific cost of modulating the area per lipid (see note (18) in Ref. [2T]). In the present 
work, we show that the terms in a, ft and C cannot be neglected with respect to others. We focus on the 
influence of a, for which we provide a physical interpretation, and we will set ft = £ = in the body of 
this article in order to simplify our discussion and to avoid adding unknown parameters. However, the 
derivation of the membrane effective Hamiltonian is presented in Sees. 11.1111.21 of our Methods part, in 
the general case where a, ft and ( are all included. 

The effective Hamiltonian of a bilayer membrane patch with projected area A p at chemical potential 
/1 can be derived from Eq. [1] For this, the effective Hamiltonians per unit projected area of the two 
monolayers are summed, taking into account the constraint that there is no space between the two 
monolayers of the bilayer, and assuming that the hydrophobic chains of the lipids are incompressible. 
This derivation is carried out in Sec. 11.11 of our Methods part. It results in an effective Hamiltonian 
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of the bilaycr membrane that depends on three variables: the average shape h = (h + + h~)/2 of the 
bilayer, the sum u of the excess hydrophobic thicknesses of the two monolayers, each being measured 
along the normal to the monolayer hydrophilic- hydrophobic interface (sec Fig. Q] and Eqs. l2"oTl2T)|) . and the 
difference 5 between the monolayer excess hydrophobic thicknesses. (The excess hydrophobic thickness 
of a monolayer is defined as the hydrophobic thicknesses of this monolayer minus its equilibrium value.) 

In the present work, we are not interested in the degree of freedom 6, which is not excited in the 
equilibrium shape of a membrane containing up-down symmetric mismatched proteins (see see Fig.[Tj3). 
Hence, in Sec. 11.21 of our Methods part, we integrate 5 out, which amounts to minimizing / with respect 
to S since our theory is Gaussian. The resulting effective Hamiltonian, which involves h and u, is given 
by Eq. [32] in Sec. 11.21 of our Methods part. In this effective Hamiltonian, the variables h and u are 
decoupled, and the part depending on h corresponds to the Hclfrich Hamiltonian [3]. Hence, our model 
gives back the Helfrich Hamiltonian if the state of the membrane is described only by its average shape 
h (sec Methods, Sec. H~3j) . 

Here, we focus on variations of the membrane thickness, i.e., on the variable u. We thus restrict 
ourselves to the case where the average shape h of the membrane is flat (see Fig. [TJ3) - in this case, we 
obtain, from Eq. [32) 



+ At V 2 u + A 2 V ■ ( uVu) H — det(didju) . (2) 

In the case where /? = £ = 0, on which the body of this article focuses, the various constants introduced 
in Eq. [5] read: 



2fi 



(3) 



K = -^(c -^ ) + k' a + j, (4) 
do 4 

K = f , (5) 

A 2 = ^(c -c E ). (7) 

In these equations, do denotes the equilibrium hydrophobic thickness of the bilayer membrane, a plays 
the part of an externally applied tension (see Methods, Sec. [5]), K a is the compressibility modulus of the 
membrane, k is its Gaussian bending rigidity, kq is the bending rigidity of a symmetric membrane such 
that 5 = 0, Co is the spontaneous (total) curvature of a monolayer, and c' is the modification of this 
spontaneous curvature due to area variations. In addition, we have introduced k' a = 4aSo/^o> which has 
the dimension of a surface tension, like K a . Note that the last three terms in Eq. [2] are boundary terms. 

In Sec. 11.21 of our Methods part, the expressions of K' a , K", A\ and A 2 are provided in the more 
general case where /3 and C are included. 

We wish to describe a membrane with an equilibrium state that corresponds to a homogeneous 
thickness. A linear stability analysis (presented in Sec. 11.41 of our Methods part) shows that the flat 
shape is stable if K a > 0, K" > 0, and 



K > -2 . (8) 
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Discussion 



Comparison with existing models 

Our model Eq. [5] has a form similar to that of the models developed in Rcfs. [T2TH5] . However, it differs 
from these previous models on several points. First, our definition of u is slightly different. Second, we 
have included the effect of an applied tension a. Finally, the various constants in Eq. [2] have different 
interpretations, and thus different values, from the ones in the existing models. Let us discuss these 
points in more detail. 



On the definition of u 

In the present work, the variable u, which is the relevant one to study membrane thickness deformations, 
is defined as the sum of the excess hydrophobic thicknesses of the two monolayers, each being measured 
along the normal to the monolayer hydrophilic-hydrophobic interface (see Eqs. [2611291 in the Methods 
section). This definition of u has the advantage of being independent of deformations of the average 
membrane shape h. 

The excess thickness variable used in Refs. [7l IT2lfT5] rT8" ] l2"2" ] [23] reads in our notations: 

h+ - h" - d Q 

u= ^ ■ (9) 

Using Eqs. [9landl25[ and working to second order, we obtain 



o - d o 
lu = u 



2 , (V«) 2 



{Vhf + 



(10) 



which shows that there is a second-order difference between 2u and our variable u. Consequently, the 
difference between the definition used in the previous works and ours regards only the term linear in 
u, i.e., the tension term, which was not included in these works. At zero applied tension, the two 
definitions are equivalent, i.e., it is equivalent to use u or 2u. Our definition of u is the right one 
for rigorously taking tension into account, because it is independent of deformations of the average 
membrane shape h: the energy stored in the variable u only comes from thickness variations. (The 
variable u of Refs. [7l[T2T[T5, 18, 22, 23] corresponds to the difference between the bilayer hydrophobic 
thickness projected along z and the non-projected equilibrium hydrophobic bilayer thickness (see Eq. [9j, 
so it is not independent of h. The second-order difference between 2 u and u, which is shown in Eq. ITO"! 
arises from this difference in projection between actual thicknesses and equilibrium thicknesses within 
the definition of u.) 



On tension 

First of all, existing models [7] [T2T[T5l [T8 l [22 ] were constructed at zero applied tension, which means a = 
in Eq. [2] To our knowledge, our work is the first where the coefficient of the term linear in u is explicitly 
related to the applied tension (see Methods, Sec. [5]) and to the tension of the Helfrich model (see Methods, 
Sec.OJ). 

In Ref. |18j . the effect of applied tension is taken into account, in so far as it changes the equilibrium 
membrane thickness of a homogeneous membrane, but without being fully implemented in the elastic 
model. Our more complete description gives back this effect on membrane thickness, when it is applied 
to the particular case of a homogeneous membrane (see Methods, Sec. [2]). 
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On the constant K' a 

In our model, the constant K' a features three contributions with different origins (see Eq. U). 

The first contribution arises from the spontaneous curvature of a monolayer and from its variation 
with the area per lipid. More precisely, the term 

£-(co - c^Sq) uV 2 u = ^-(c - c^So) [V • (uVu) - (Vu) 2 ] (11) 

appears when one constructs the membrane model starting from a monolayer Hamiltonian density such 
as Eq. [1] This term was first introduced in Ref. [12], and it was then included in Refs. [T3lll4j . 

The second contribution, k' a , arises from a, i.e., from the term in (VE) 2 introduced in Eq. [TJ This 
term was not included in Refs. [12H14) . which started from a second-order expansion of the effective 
Hamiltonian per lipid molecule involving only the curvature and the area per lipid. However, a gradient 
of area per lipid (or, equivalently, of the thickness) in a monolayer has an energetic cost of its own. Indeed, 
a greater part of the hydrophobic chains is in contact with water when a thickness gradient is present 
(see Fig. [2]). The associated energetic cost is given by the interfacial tension 7 of the hydrocarbon- water 
interface, which is of order 40-50 mN/m [241,25!. Such a term is often accounted for in microscopic 
membrane models (see, e.g., Ref. [26]). In the case of a symmetric membrane (u + = u~ = u/2) with flat 
average shape, the surface of the hydrocarbon- water interface is increased by a factor [1 + (Vu) 2 /8] for 
each monolayer (see Fig. [2]). Thus, to second order, the associated energetic cost per unit projected area 
is 7(Vu) 2 /4. Note that other physical effects, e.g., the elasticity of the chains, may yield contributions to 
the term in (VS) 2 . However, if we restrict to the simple term arising from interfacial tension, we obtain 

fc; = |«25mN/m. (12) 

Finally, the third contribution, cr/4, arises from the (macroscopic) externally applied tension. The 
tension of a vesicle can rise only up to a few mN/m before it bursts (see, e.g., Ref. [H]). Hence, according 
to our estimate of k' a in Eq. [TU we expect er/4 <C k' a . 

In the seminal article Ref. [7], where the membrane model was constructed by analogy with liquid 
crystals, a term in (Vu) 2 , interpreted as arising from tension, was included in the effective Hamiltonian. 
However, its effect was neglected on the grounds that the value of its prefactor made it negligible with 
respect to the other terms. The value of this prefactor was taken to be that of the tension of a monolayer 
on the surface of a Plateau border [27]. The model introduced in Ref. [7] was further developed and 
analyzed in Refs. [i"8]l22]. where the same argument was used to neglect the term in (Vw) 2 . 

However, our construction of the membrane effective Hamiltonian shows that the microscopic tension 
involved through k' a arises from local variations in the area per lipid. This stands in contrast with the 
case of the Plateau border, where whole molecules can move along the surface and exchange with the 
bulk, yielding a smaller value of the tension. Ref. [27] stresses that the measured tension of a Plateau 
border is valid for long-wavelength fluctuations, but that it is largely underestimated for short- wavelength 
fluctuations (less than 10 nm) which involve significant changes in area per molecule. 

Including the tension of the hydrocarbon-water interface instead of that of the Plateau border is a 
significant change, given that the former is of order 40-50 mN/m [2~4l,l2"5] . while the latter is of order 
1.5-3 mN/m [7l H8l[22l[27] . In Refs. (T8II22"] . it is shown that the effect of the term in (Vu) 2 is negligible if 

« , (13) 

«o 

where we have used our own notations of the prcfactors of the terms in (Vu) 2 , u 2 and (V 2 u) 2 . In the case 
of DOPC, taking K'l = k/4 and using the values of the membrane constants [28], this condition becomes 
K' a <C 28 mN/m. While this is well verified if K' a corresponds to the tension of the Plateau border, it is 
no longer valid within our model. 
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Our model is the first that includes all contributions to K' a , in particular the one arising from interfacial 
tension. Besides, in Sec. ll.2l of our Methods part, we show that /3 is also involved in K' a , which emphasizes 
the complexity of constructing a continuum model to describe membrane elasticity at the nanoscale: many 
terms involved in the expansion of the effective Hamiltonian cannot be neglected a priori. 

In the following, we will analyze numerical and experimental data, looking for evidence for the presence 
of k' a , and comparing the relative weight of the different contributions to K' a . 

On the value of K' c [ 

We have obtained K'^ = kq/A (see Eq. [3]) , where «o is the bending rigidity of a symmetric membrane 
such that (5 = 0. The clastic constant kq is related to the bending rigidity re of the Helfrich model (see 
Methods, Sec. [OJ through 

K = KQ -^(co-c^ ) 2 . (14) 

The difference between «o and re arises from integrating out S (see Methods, Sec. II. 2p . In the previous 
models, this procedure was not carried out, as one focused directly on the symmetric case 6 = 0. All 
previous models thus made the approximation K'J = re/4 [7] IT2T4l"4^ [l"8 l l2"2"] . 

In addition, in Sec. 11.21 of our Methods part, we show that £ is also involved in K", which stresses 
further the possible importance of such terms in order to describe membrane elasticity at the nanoscale. 

On boundary terms 

The boundary terms correspond to the last three terms in Eq. [2] When one wishes to describe the local 
membrane deformation due to a transmembrane protein, boundary terms play an important part, as their 
integral on the contour of the protein contributes to the deformation energy. The first two boundary 
terms are the same as in Refs. [T2H14] . However, even at vanishing applied tension, we have K' a ^ —2A2, 
contrary to the previous models [H], due to the presence of k' a . We have also accounted for the Gaussian 
bending rigidity re, as in Ref. [15] : it yields the third boundary term. 

Again, the situation is more complex when j3 is included, as the expressions of A\ and A2 then feature 
extra terms linear in f3 (see Eq. |37]in Sec. 11.21 of our Methods part). 

On lipid tilt 

Several membrane models including lipid tilt in addition to average shape deformations and/or thickness 
deformations have been elaborated [21.23, 26 , 29 3Tj. These models provide improvements with respect 
to the Helfrich model, yielding better agreement with numerical data on bulk membranes [23 ] l31 j . 

Our model docs not include lipid tilt because we focus on local thickness deformations, and especially 
on comparison to experimental and numerical data regarding deformations induced by mismatched pro- 
teins. While it would be interesting to include this extra degree of freedom, it would imply introducing 
several membrane parameters, which would make comparison to mismatch data impractical. 

Not taking tilt into account means that we are effectively integrating out this degree of freedom 
through coarse-graining. More precisely, the elastic coefficients of a more detailed membrane model, which 
would include tilt as an extra degree of freedom, would be rcnormalized by integrating out tilt. This means 
that tilt is included within the elastic coefficients of our membrane model. In addition, the interaction 
energy between the membrane and a mismatched inclusion (see, e.g., Eq. ITSl) . and, consequently, the 
effective boundary conditions at the inclusion boundary, may involve tilt (see, e.g., Ref. [21]). In this 
interaction energy, tilt can be integrated out in the same way as in the bulk membrane energy. Hence, 
we are not losing any part of the elastic energy by disregarding the tilt degree of freedom. However, it is 
not impossible that a model including tilt truncated at second order could prove more efficient (e.g., have 
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a wider domain of validity at short wavelengths) than one truncated at the same order and disregarding 
tilt. 

Comparison with numerical results 

As numerical simulations become more and more realistic, they start providing insight into the behavior 
of systems on the microscopic scale where direct experimental observation is difficult. Lipid membranes 
(with or without inclusions) are no exception. Over the last decade, several groups have simulated bilayer 
systems over length- and time-scales long enough to give access to the material constants relevant for 
nanoscale deformations. These simulations provide interesting tests for theoretical models describing 
membrane elasticity at the nanoscale. We will compare the predictions of our model to recent numerical 
results in this Section. All the numerical results we will discuss have been obtained at zero applied 
tension. Hence, throughout this section, we take a = 0. This implies that our definition of the membrane 
thickness is equivalent to that considered in the original numerical works (see the discussion above on 
the definition of u). 

Fluctuation spectra 

Using numerical simulations, one can measure precisely the fluctuation spectra of the average height and 
the thickness of a bilayer membrane [14, 16 ,32,33 . Microscopic protrusion modes, occurring at the scale 
of a lipid molecule, contribute to these spectra. While they are not described by continuum theories, it is 
possible to consider that they are decoupled from the larger-scale modes fl^lTS] . By fitting the numerical 
spectra to theoretical formulas, it is possible to extract the numerical values of the membrane constants 
involved in the continuum theory. In our framework, the fluctuation spectra of the average height of 
the membrane give access to the Helfrich bending rigidity k, while those regarding the thickness of the 
membrane give access to K a , K' a and K^. 

We have reanalyzed the height and thickness spectra presented in Rcfs. [16,32, 33] using the fitting 
formulas in Refs. [T4fl6j (see Eq. 32 of Ref. [14]) and the method described in Ref. [14], except that we did 
not assume that K 1 ^ = k/4, in order to include the possible effect of the difference between k and kq (see 
Eg. |3"3"]) . and of ( (see Eg. [33]) . Our results were similar to those obtained in Refs. [TlirTo] assuming that 
A'" = k/4, and we obtained no systematic significant difference between K'^ and k/4, which means that 
the corrections to A"" predicted by our model are negligible in these simulations. This gives a justification 
for focusing only on the correction to K' a , as we do in this article. Besides, we obtained K' a < from all 
the fits, as reported in Refs. [T4" l ll6j. and we checked that all the values obtained for K' a comply with the 
stability condition Eq. [5] 

Deformation profiles close to a mismatched protein 

In Refs. [14f(16j. the thickness profile of a membrane containing one cylindrical inclusion with a hydropho- 
bic mismatch has been obtained from coarse-grained numerical simulations. Comparing the average nu- 
merical thickness profiles to the equilibrium profiles predicted from theory is a good test for our model, 
in particular to find clues for the presence of k' a . 

Let us denote the radius of the protein by ro, and its hydrophobic length by b. the mismatch originates 
from the difference between I and the equilibrium hydrophobic thickness do of the membrane. The 
equilibrium shape of the membrane, which minimizes its deformation energy, is solution to the Eulcr- 
Lagrange equation associated with the effective Hamiltonian density in Eq. O We write down this 
equilibrium shape explicitly in Sec. !3.1l of our Methods part. In order to determine it fully, it is necessary 
to impose boundary conditions at the edge of the inclusion, i.e., in r = ro- There is a consensus on the 
assumption of strong hydrophobic coupling u(ro) — u$ = £— do, as it costs more energy to expose part of 
the hydrophobic chains to water than to deform the membrane, for typical mismatches of a few A. Note 
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that, with our definition of u, the condition u(ro) = uq = £ — do is valid to first order, while it is exactly 
valid with the definition of Refs. [71 [TW51[T51[^1^5] (see Eqs. H [TO]). This difference arises from the fact 
that our u is not projected along z (sec Fig. |TJ), which makes it fully independent of h. Given that the 
clastic energy is known to second order, the equilibrium membrane shape resulting from its minimization 
is known to first order, so it is sufficient to use boundary conditions to first order. Hence, such differences 
are not relevant for the present study and will not be mentioned any longer. 

However, there is some debate about the second boundary condition in rg (see, e.g., Ref. |14j). which 
regards the slope of the membrane thickness profile. Traditionally, one either assumes that the protein 
locally imposes a fixed slope to the membrane |18ll22j . or minimizes the effective Hamiltonian in the 
absence of any additional constraint, which amounts to considering that the system is free to adjust its 
slope in r = r |12H16j . In Sec. 13. II of our Methods part, we present the equilibrium profiles for these two 
types of boundary conditions. The actual boundary condition depends on the interactions between the 
protein and the membrane. In a quadratic approximation, these interactions gencrically give rise to an 
effective potential f s favoring a slope sq in tq: 



where k s is an effective rigidity, while u' denotes the derivative of the membrane thickness profile u with 
respect to the radial coordinate r. Two a priori unknown parameters, A; s and Sq, are associated with this 
effective potential. The "free-slope" boundary condition (also called "natural" boundary condition [12[ 
[14]) is recovered in the limit k s — > 0, which is appropriate if f s is negligible with respect to the energetic 
contributions in /. Conversely, if k s — > oo, the protein locally imposes the fixed slope sq. If the interactions 
between the protein and the membrane lipids are sufficiently short-ranged, the protein cannot effectively 
impose or favor a slope on the coarse-grained membrane thickness profile. For instance, in the numerical 
simulations of Rcfs. |14H16) . the interactions between the protein and the membrane lipids are of similar 
nature and of similar range as those between membrane lipids. Thus, we will choose the free-slope 
boundary condition in our analysis of this data. This choice was already made in Refs. |14H16j . A 
practical advantage of this boundary condition is that it does not introduce any unknown parameter in 
the description. 

The membrane model of Refs. [T4HT6] is very similar to ours, except that k' a = 0. It was shown in 
Ref. |16j that this model can reproduce very well the numerical results, provided that the spontaneous 
curvature is adjusted for each deformation profile (see Fig. [3]). In Ref. |16j . the adjusted "renormalized 
spontaneous curvature", denoted by Co, was found to depend linearly on the hydrophobic mismatch 
Mo |16| , as shown in Fig. |4] In our model, the equilibrium profile corresponding to the free-slope boundary 
conditions (see Eqs. I46land l53]) involves k' a . We show in Sec. 13.11 of our Methods part that the quantity 



then plays the part of the renormalized spontaneous curvature of Ref. |16| in the equilibrium profile. This 
quantity is linear in uq: our model, and more precisely the presence of a nonvanishing k' a , thus provides 
an appealing explanation for the linear dependence observed in Ref. j!6j . 

Using a linear fit of the data of Ref. [16] (see Fig. [4]), together with Eq. [16] and the value k = 
2.3 x 10~ 20 J extracted from the spectra in Ref. [15] , we obtain k' a = 13mN/m. 

It is interesting to compare this value to K' a , which is obtained from the fluctuation spectra in 
Ref. [IB]: K' a = — 9.2mN/m. This shows that the contribution of k' a to K' a is important. Besides, we 
may now estimate the contribution to K' a that stems from the monolayer spontaneous curvature (see 
Eq. [4]): — k (cq — c' T l0 )/do = K' a — k' a = — 22mN/m. Using values from the fluctuation spectra in 
Ref. [16], this yields £ « — 6 A for the algebraic distance from the neutral surface of a monolayer to the 
hydrophilic-hydrophobic interface of this monolayer (see Methods, Sec. 0] for the relation between £ and 



f s = k s (u'(r ) - so) 



(15) 



co = c H Mo , 



(16) 
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In Rcf . |15j , a different coarse-grained molecular simulation model was used to obtain the equilibrium 
membrane thickness profiles for cylindrical inclusions with two different hydrophobic thicknesses, one 
yielding a positive mismatch and the other a negative one, and with seven different radii ro- These 
profiles are presented in Figs. 6 and 7 of Ref. [TS], except those corresponding to the inclusions with 
largest radii (5.25 nm), but this data was communicated to us by one of the authors of Ref. [15]. We 
fitted each of these numerical profiles to the analytical equilibrium profile Eq.06]with prcfactors ^4±(0, Co) 
(see Eq. [54]). using cq as our only fitting parameter, in the spirit of Rcf. [16]. We found that cq does not 
depend on the radius of the inclusion, but that it depends significantly on the mismatch (sec Fig. EK). 
This is in good agreement with the predictions of our model (see Eq. ITS]) . For each of the two values of 
wo, we have averaged cq over the seven results corresponding to the different inclusion radii. The line 
joining these two average values of cq as a function of Uq is plotted in Fig. [SJ3. Using Eq. [JJ)] and the 
value k = 1.4 x 10~ 19 J [HUTS], the slope of this line yields k' a = 36mN/m: the order of magnitude of 
this value is the same as the one obtained from the data of Ref. PJj] . 

Again, we can compare this value to K' a , which is obtained from the fluctuation spectra in Refs. [T41I15) : 
K' a = — 11.9mN/m. Hence, the contribution of k' a to K' a is important here too. We also obtain — kq(cq — 
c S )/d = K' a -k' a = -48mN/m, and £ w -3 A. 

In Ref. [15] , the shortcomings of the model that disregards k' a are explained by the local variation of 
the volume per lipid close to the protein. It is shown in Ref. |15j that including this effect yields 

c = c - —v{r ) , (17) 
vo 

where vq is the bulk equilibrium volume per lipid, while vq +v(ro) denotes the volume per lipid in r = tq. 
However, the predicted linear dependence of Co in v(r )/vQ is not obvious: in Fig. [S] we rather see two 
groups of points (one for each value of Uq) than a linear law. In other words, the data of Ref. [15] is 
more consistent with a value of Co that depends only on % and not on v (or r ), in agreement with the 
predictions of our model (see Eq. [To]) . In Ref. PJj], local modifications of the volume per lipid close to 
the inclusion were investigated too, as well as local modifications of the nematic order, of the shielding 
of the hydrophobic membrane interior from the solvent, and of the overlap between the two monolayers. 
None of these effects was found to explain satisfactorily the linear dependence of Co versus uo |16| . 

To sum up, our model can explain the dependence of Co in uo observed in the numerical results 
of Refs. [ljjinji] as a consequence of the presence of k' a . Our explanation does not involve any local 
modification of the membrane properties, in contrast with those proposed in Refs. |15U16j . Furthermore, 
the order of magnitude we obtain for k' a from the data of Refs. |15U16| is in agreement with our estimate 
in Eq. [JJ 

Comparison with experimental results 

The antimicrobial linear pentadecapeptide gramicidin (see [8] for a review) is a very convenient exper- 
imental system to probe membrane elasticity on the nanoscale. In lipid membranes, two gramicidin 
monomers (one in each monolayer) associate via the N-terminus to form a dimeric channel, stabilized 
by six intermolecular hydrogen bonds. The channel being large enough for the passage of monovalent 
cations, conductivity measurements [9] can detect its formation and lifetime, which are directly influ- 
enced by the membrane properties. Indeed, while the monomers do not deform the membrane, the 
dimeric channel presents a hydrophobic mismatch with the membrane, so that dimer formation involves 
a local deformation of the bilayer. The gramicidin channel can therefore act as a local probe for the 
bilayer elasticity. Furthermore, the gramicidin channel can be considered as up-down symmetric and 
cylinder-shaped, which makes it convenient for theoretical investigations. 

Data on gramicidin channels originally motivated theoretical investigations on membrane models 
describing local thickness deformations [7] ITTlfT3"] . Such data now provides a great opportunity to test 
any refinement of these models. We will compare our model to the data of Ref. [17] regarding the lifetime 
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of the gramicidin channel as a function of bilayer thickness, and then to the data of Ref. [H] regarding 
the formation rate of the gramicidin channel as a function of bilayer tension. 

In order to compare the predictions of our model to experimental data regarding the gramicidin chan- 
nel, it is necessary to make some assumptions about the boundary conditions at the edge of the channel, 
i.e., in r — ro- As discussed in the previous section, we will assume strong hydrophobic coupling, i.e., 
u(ro) = Uo = £ — do, but determining the boundary condition on the slope of the membrane thickness 
profile is trickier as it depends on the interactions between gramicidin and the membrane lipids. In 
previous analyses |18[I34]. the fixed-slope boundary condition was favored as giving the best agreement 
with experimental data. However, different values of the fixed slope were obtained in these studies. In 
addition, recent all-atom simulations of gramicidin channels in lipid bilayers indicate that the membrane 
thickness profile is complex in the first lipid shell around the channel, due to specific interactions, and 
that beyond this first shell, no common slope exists for the different membranes investigated |35j . Given 
the difficulty to determine the actual effective boundary condition associated with the slope of the mem- 
brane thickness profile, we will adopt the free-slope boundary condition, which has the advantage not to 
introduce any unknown parameter in the analysis, but we will also compare our results to those obtained 
with the more traditional fixed-slope boundary condition. 

Analysis of the experimental data of Elliott et al. |17] 

It was shown in Ref. [22] that the detailed elastic membrane model introduced in Ref. [7] yields an 
effective linear spring model as far as the membrane deformation due to gramicidin is concerned [2"2"ll34j : 
the energy variation F associated with the deformation can be expressed as F = Huq, where H is the 
effective spring constant, while u$ is the thickness mismatch between the gramicidin channel and the 
membrane. This linear spring model was validated by comparison with experimental data on the lifetime 
of the gramicidin channel, measured as a function of bilayer thickness ( |17[l36j. summarized in [34]) and 
function of the channel length [57] , 

We will here focus on the data concerning virtually solvent-free bilayers, i.e., membranes formed using 
squalene. The elasticity of membranes containing hydrocarbons should be different: for instance, a local 
thickness change of the membrane could be associated with a redistribution of the hydrocarbons. (In 
this, our analysis differs from that of Ref. [14], where all the data of Ref. [17] was considered. Another 
important difference with the analysis conducted in that reference is that we use experimental values of 
the membrane parameters, which are quite different from the values coming from numerical simulations.) 
In Ref. [34], the effective spring constant H of the membrane was estimated from data of Ref. [17] on 
gramicidin channel lifetime for three bilayers formed in squalene with monoglycerids that differed only 
through their chain lengths: the different thicknesses of these membranes yield different hydrophobic 
mismatches with a given type of gramicidin channels. The value H cxp = 115 ± 10 mN.m -1 was obtained. 

In Sec. l3.2l of our Methods part, we use our model to calculate the deformation energy of the membrane 
due to the presence of a mismatched protein. Both in the case of the free-slope boundary condition, and 
in the case where the gramicidin channel locally imposes a vanishing slope, this deformation energy can 
be expressed as a quadratic function of the mismatch uq. The prefactor of Uq in the deformation energy 
F corresponds to the effective spring constant of the system. Thus, although our model is different from 
the one of Refs. [7JIS1122], it also yields an effective linear spring model. This is not surprising since we 
are dealing with the small deformations of an elastic system. However, the detailed expressions of our 
spring constants as a function of the membrane parameters (see Eqs. 1591 and I65|) arc different from those 
obtained using the model of Refs. [2IIB1122], due to the differences between the underlying membrane 
models. In particular, in our model, k' a is involved in H, through K' a . Our aim will be to find out which 
value of k' a gives the best agreement with the experimental value of H. 

Using Eqs. 2] [5] and and neglecting the difference between k and Ko, Eqs. [55] and [55] show that 
H depends on the elastic constants ft, R and Co involved in the Helfrich model, on K a , on c £oj which 
corresponds to the spontaneous curvature variation with the area per lipid, on do, on the radius ro of 
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the gramicidin channel, and on k' a . There is, to our knowledge, no direct experimental measurement of 
CgEo available, but, as shown in Sec. 2] our Methods part, we have CqSo = A' q £/k, where £ denotes the 
algebraic distance from the neutral surface of a monolayer to the hydrophilic-hydrophobic interface of 
this monolayer (see Eg. 1731 neglecting the difference between k and kq). Hence, in order to calculate the 
spring constant, we need values for k, R, Co, K a and £, in the precise case of monoolein membranes. 

In Rcf. [38] . the elastic constants k, R and cq were measured in a monoolein cubic mesophase, both at 
25°C and at 35°C. The positions of the neutral surface and of the hydrophilic-hydrophobic interface were 
estimated on the same system in Ref. [39], but these results were flawed by a mathematical issue, which 
was corrected in Ref. [50]. This correction yielded other corrections on cq, and on the ratio R/k |41) . 
These results regard a cubic phase, where the membrane is highly deformed with respect to a flat bilayer: 
the values of the various constants should be affected by the strains present in this phase. In another 
work [42] . the constants of monoolein are determined in a highly hydratcd doped Hn phase, where the 
strains are better relaxed. However, these measurements were carried out at 37°C, while the experiments 
of Rcf. [T7] that we wish to analyze were performed at 23°C. Given that the data of Refs. [55II3T?] include 
the most appropriate temperature, while the ones of Ref. [42) correspond to the most appropriate phase, 
we will present results corresponding to both sets of parameters. Finally, the experimental value of K a 
for monoolein is provided by Ref. [27]. 

In Table Q] we present the results obtained for the spring constant H of monoolein bilayers, using 
the different experimental estimates of the membrane constants. The main difference between parameter 
sets 1 and 2 is the value and the sign of R [351111] ■ However, R is involved in H only in the free-slope 
case (see Eqs. [591 and [65]) : the 3% difference between the values of H obtained with parameter sets 1 
and 2 stems only from the difference on Co, while the 12% difference between Hf obtained with data sets 
1 and 2 contains an important contribution from R. The constants in parameter set 3, corresponding to 
Ref. [52], are significantly different from those of Refs. [3"51l41j . which yields a 30% difference on H and a 
20% difference on Hf. We also note that, as the value of the algebraic distance from the neutral surface 
to the hydrophilic-hydrophobic interface of a monolayer is very small compared to the other length scales 
involved (£ = —0.3 A [40]), the contribution of CqEo to H is negligible (it is of order 1%). 

Let us now discuss the results given by our model, in the case of the free-slope boundary condition (see 
Tabled]). The spring constants Hf obtained assuming that k' a = are about three times smaller than the 
experimental value H oxp = 115 ± 10 mN.m -1 (see line 1 of Table [TJ. (This result is very similar to that in 
Ref. [34] , which illustrates that accounting for monolayer spontaneous curvature and for boundary terms 
does not change much the value of H.) However, Hf reaches the experimental value for k' a « 25mN/m 
for all three parameter sets (see line 2 of Table [lj. Hence, for free-slope boundary conditions, the presence 
of k' a , with an order of magnitude consistent with Eq. I12[ improves the agreement between theory and 
experiment. 

We may compare these values of k' a to the contribution to K' a that originates from the monolayer 
spontaneous curvature (see Eq. [4j: — k (c — CqSo)/^- We estimate the value of this contribution to be 
between 0.26 and 1.2mN/m, depending on which set of parameters is chosen. This is positive and much 
smaller in absolute value than the estimates obtained from the numerical data of Ref. [16] and of Ref. [15] : 
here, the neutral surface of a monolayer and its hydrophilic-hydrophobic interface are very close, while £ 
seemed to be significant in the numerical simulations. In addition, the contribution of membrane tension 
to K' a , namely, cr/4, cannot exceed about 1 mN/m. In the case of the free-slope boundary condition, our 
results imply that k' a should be the dominant contribution to K' a for the membranes studied in Ref. [17] . 

Let us now discuss the results obtained for the zero-slope boundary condition, which was investigated 
in Ref. [34]. For the zero-slope boundary condition, the values obtained for Hq assuming that k' a = 
are in quite good agreement with the experimental value H exp = 115 ± 10 mN.m -1 obtained in Ref. [34j 
from the data of Ref. [17], for all the data sets we used (see line 3 of Table [T]): hence, k' a seems negligible 
if zero-slope boundary conditions are assumed. However, there is no justification to assume that the 
gramicidin channel locally imposes a vanishing slope. 
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Analysis of the experimental data of Goulian et al. |18| 

While the experiments cited in the previous Section dealt with discrete changes of the hydrophobic mis- 
match obtained by varying membrane composition, Goulian et al. j!8) measured the gramicidin channel 
formation rate / in lipid vesicles as a function of the tension a applied through a micropipette. As the 
tension is an externally controlled parameter that can be changed continuously for the same gramicidin- 
containing membrane, this approach can yield more information, and it has the advantage of limiting 
the experimental artifacts associated to new preparations. To date, the experiment in Ref. |18) remains 
the most significant in the field and should serve as a testing ground for any theoretical model. We will 
therefore discuss in detail the data and its interpretation by the original authors [18l[22j as well as in 
terms of our model (see Eq. [2]). 

Within experimental precision, the data of Ref. |18j can be described by a quadratic dependence: 

In/ = 5 (a) = C + C 1( j + C 2 a 2 . (18) 

Given that In / is a linear function of the energy barrier associated with the formation of the gramicidin 
dimcr, it is a sum of a chemical contribution, including, e.g., the energy involved in hydrogen bond 
formation, and of a contribution arising from membrane deformation due to the dimer (monomers do 
not deform the membrane) [18J . The latter contribution arises from the hydrophobic mismatch between 
the membrane and the dimer, and it depends on the applied tension a, since the membrane hydrophobic 
thickness depends on a (see Eq. 03] in Sec. [2] of our Methods part). Expressing the deformation energy 
F of the membrane due to the presence of the dimer gives a theoretical expression for the tr-dependent 
part of In/. In our model, K' a features a contribution coming from a (see Eq. [4]). However, this term is 
negligible, given that a/4 -C K a K"/ d$ (see Eq.H5]). for realistic tension values (a few mN/m at most), 
and for the experimentally measured values of the membrane constants [28] . This enables us to disregard 
it. Then, our quadratic elastic membrane model simply gives a quadratic dependence of F on a, in 
agreement with the form of Eq. 1181 Comparing the experimental values of C\ and Ci to those predicted 
by theory provides a test for theoretical models [18]. (Note that, if the cr-dependent contribution to K' a 
is included, the expression of the er-dependent part of In / is no longer simply quadratic in a. However, 
we explicitly verified that including this contribution yields a negligible change to the relation between a 
and In /, for realistic values of the parameters). 

Since the coefficients C\ and Ci arise from membrane elasticity, they are common to all the vesicles 
studied in Ref. |18j . which have the same lipid composition. Conversely, the baseline Co depends on 
parameters such as the concentration of gramicidin molecules, so it can take a different value for each 
of the twelve vesicles studied in Ref. [T5] . A global fit to the data of Ref. [TH] using Eq. [15] involves 
minimizing the goodness-of-fit function 

x 2 = E( ln ^--^)) 2 ' ( 19 ) 

3 

where the index j runs over all the experimental points, with fitting parameters Ci, Ci, Cg , k = 1, ... ,12. 
The baseline Cq is then subtracted from each of the twelve curves. All the data is plotted in the 
same graph in Fig. [7J The best global fit, corresponding to C\ = 0.74 ± 0.07 (mN/m) ~ and C 2 = 
-0.090 ± 0.015 (mN/m) , is shown on Fig. [7] as the dotted (black) line. (It should be noted that 
the values obtained by fitting the individual curves are much more scattered: C\ ranges from 0.4 to 
1.5 (mN/m) -1 and C 2 from -0.3 to (mN/m) -2 .) 

In Ref. [18] , the authors used published values of the material constants to calculate C\ and Ci in the 
framework of their elastic model [35], based on that of Ref. [7J. Using fixed-slope boundary conditions, 
they reported good agreement with the experimental data for a reasonable value of the unknown slope s 
(s = 0.3). However, we need to raise the following points: 
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1. There was a mistake in their implementation of the formula of Ref. |22] giving C\ and C2 as a 
function of the material constants. More precisely, we found that a factor of 2 was missing in the 
expression of C\ and a factor of 4 was missing in that of C2 in the implementation of the formula 
of Ref. [35]. This was confirmed by Mark Goulian (private communication). The actual values of 
C\ and C2 obtained using the same values of the constants as in Ref. [18] are in fact quite far from 
those corresponding to the best fit of the experimental data, as shown by the dashed green line in 
Fig. [Jj (sec also Fig. M and Table [2]) . 

2. The estimates for the elastic constants used in Ref. [TB] are somewhat different from more recent 
and more widely accepted values. Henceforth, we will use the following parameters, for a DOPC 
membrane: d = 2.7nm [18], K a = 265mN/m, k = 8.5x 1CT 20 J [28], c = -0.132nm- 1 [43], and 
the dimensions of a gramicidin channel: ro = 1 nm, £' = I + S = 2.3 nm [18]. Implementing these 
more recent values in the model of Ref. [22j does not yield a better agreement with experiment, as 
shown by the dashed-dottcd (blue) line in Figure [7] (see also Fig. [8] and Table [2]). 

A somewhat better agreement with the experimental data is obtained when taking s = instead of 
s = 0.3 for the fixed slope (see Figs. [7] and [8] and Table [2]). However, the downward inflection of the 
experimental curves at high a is not adequately described for any value of ,s. In fact, C2 is independent 
of s, and its absolute value given by the elastic model is 15 times smaller than the experimental one (see 
Table [5]). We conclude that the elastic model of Refs. [7JH2] does not satisfactorily describe the data of 
Ref. |18j regarding the lifetime of the gramicidin channel under tension. 

In Sec. 13.21 of our Methods part, we calculate the deformation energy F in the framework of our 
model, both for the fixed-slope boundary condition and for the free-slope boundary condition. The 
resulting expressions of C\ and C2 are given by Eqs. [61] [62] [67] and [68] In order to see which values of 
k' a and which boundary conditions give the best agreement with the experiments of Ref. |18| , we present 
a plot of the goodness-of-fit function \ 2 ( see Eq. H2) in a (Cii C2) graph in Fig. [5] On this graph, we 
have plotted the trajectories obtained from our model in the (Ci, C2) plane when varying k' a , for s = 0, 
for ,s = 0.3 (as in Ref. [IB]), and for the free-slope boundary condition. 

In order to obtain numerical values of C\ and C2 from Eqs. [61] [62] [67] and [68] we used the above- 
mentioned parameter values, and the estimate R = —0.8k [16] . Finally, we estimated CqSo through the 
relation CqSo = K<i£,/k (see Eq. [73] in Sec. 0] of our Methods part). For this, the algebraic distance 
£ from the neutral surface of a monolayer to the hydrophilic-hydrophobic interface of this monolayer 
was estimated by first determining the position of the pivot surface from the data of Ref. [43] , and by 
calculating the distance between it and the neutral surface [44]: we found £ « —0.5 A. Here again, 
the neutral surface is close to the hydrophilic-hydrophobic interface. For the sake of simplicity, we took 
Cq = 0, and we checked that the results were not significantly different when taking £ = —0.5 A. 

The ingredient in our model that can change significantly the results is k' a (Note that the values of 
C\ and C2 corresponding to k' a = arc very close to those obtained using the model of Ref. [18] with our 
values of the parameters, as shown in Table [2] This illustrates again that the influence of boundary terms 
is quantitatively small.) Fig. [8] shows that the experimental value of C\ can be explained by our model. 
In addition, the values of k' a that minimize x 2 , i.e., that give the best agreement with the experimental 
data of Ref. [T5], are between and 50mN/m, depending on the boundary condition chosen, as shown 
in Table [3] This range of values of k' a is reasonable. 

For the free-slope boundary condition, the best agreement with the experimental results is obtained 
for k' a ~ 40mN/m (see Tableland Fig. [8]). The order of magnitude is the one expected from k' a ~ 7/2. 

Let us now discuss the results obtained for the fixed-slope boundary condition, which is used in 
Ref. [T5]. For a fixed slope s = 0, the best agreement with the results of Ref. [T5] analyzed with the 
complete quadratic fit is obtained for k' a = 0. Conversely, for s = 0.3, the best agreement is obtained 
for k' a sa 40mN/m, which is similar to the result obtained the free-slope case (see Table [3] and Fig. [8]). 
Hence, in the case of the fixed-slope boundary condition, the conclusions depend a lot on the value of s 
that is chosen. 
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In all cases, the absolute values of Ci we obtain remain much smaller than the one that matches best 
the experimental results, which is Ci = —90.0 x 10~ 3 (mN/m) -2 (see Fig. [7]). This can be seen in Fig. [SJ 
as well as in Table [3] Hence, with our model, as with the one of Ref. [TH] , it seems impossible to explain 
the experimental value of Ci . Our model predicts that Ci is proportional to the effective spring constant 
H of the membrane discussed in the previous Section (see Eqs. I59l and l65[) : it is thus quite unexpected to 
have a good agreement with the experimental values of H but not with those of Ci ■ This disagreement 
on C2 could come either from a shortcoming of the model or from an undetected systematic error in 
the experimental data. The importance of Ci is largest at highest tensions, as it is C2 which gives the 
curve its concavity, and it should be noted that the maximum applied tension a is around 4.5 mN/m in 
Ref. [18], which is comparable to the rupture threshold of 3 — 10 mN/m [28]. The membrane properties 
may be affected at such high tensions in a way that is no longer well described by standard clastic models. 
It would be interesting to have more experimental data on the behavior of gramicidin channels under 
tension to see if this unexpected value of Ci persists. 

Following the hypothesis that high tensions are problematic, we performed a linear fit of the data 
of Ref. [18] (i.e., a fit with C2 = 0), keeping only the points corresponding to a < 2 mN/m: this yields 
Ci = (0.62 ± 0.05) (mN/m) -1 (sec Fig.©. In TableU we list, for different boundary conditions, the 
value of k' a which gives C\ = 0.62 (mN/m) -1 , and the value of C2 obtained from our model for this k' a . 
These values correspond to those that give the best agreement between our model and the linear fit to 
the low-tension data of Ref. [18] presented in Fig. [9] Table [4] shows that the values of k' a that yield the 
best agreement with the experimental data have a similar order of magnitude as those obtained above 
with the full quadratic fit (see Tabic [3]), remaining below 100 mN/m. Again, these values depend a lot 
on s for fixed-slope boundary-conditions. (For instance, the slope s = —0.17 is consistent with k' a = 
(see Table0]). However, there is no a priori reason for assuming that k' a = 0.) 

Again, we may compare our estimates of k' a (see Tables [3] and [4]) to the term — kq(cq — Cq£o) / do, which 
also contributes to K' a : here, —kq(cq — c Y, )/do = —0.76 mN/m. This is much smaller in absolute value 
than the corresponding estimates obtained from the numerical data of Ref. [16] and of Ref. [15]: here, as 
in the membranes studied in Ref. j!7] , the neutral surface of a monolayer and its hydrophilic- hydrophobic 
interface are very close, while £ seemed to be of a few A in the numerical simulations. We note in passing 
that this hints at a relevant difference between simulated membranes and real membranes. Besides, in the 
case of the free-slope boundary condition, our results imply that k' a should be the dominant contribution 
to K' a for the membranes studied in Ref. [18] , as for those of Ref. [17] . 

Hence, for the free-slope boundary condition, our analyses of the numerical data of Ref. |16] and of 
Ref. [TS] , and our analyses of the experimental data of Ref. [T7] and of Ref. [T5] all converge toward a 
value of a few tens of mN/m for k' a , which is of the order of magnitude expected if k' a = 7/2. Conversely, 
for the fixed-slope boundary condition, the value of k' a is coupled to that of the slope s. 

Conclusion 

We have put forward a modification of membrane elastic models used to describe thickness deformations 
at the nanoscale. We have shown that terms involving the gradient (and the Laplacian) of the area 
per lipid contribute to important terms of the effective Hamiltonian of the bilayer membrane. We have 
reanalyzed numerical and experimental data to find some signature of the presence of these terms. Using 
the free-slope boundary condition at the boundary of the mismatched protein, we have obtained consistent 
results showing that the term stemming from the gradient of the area per molecule has a prefactor k' a in 
the range 13 — 60 mN/m. Such values are consistent with the idea that this term involves a significant 
contribution of the interfacial tension 7 between water and the hydrocarbon-like hydrophobic part of the 
membrane. Indeed, this contribution should yield k' a = 7/2 « 25 mN/m. 

Interestingly, our analysis of the experimental data from Ref. |18] has shown that these nice experi- 
mental results were not as well understood as assumed in the literature. Hence, it would be interesting 
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to have more data on the behavior of gramicidin channels in membranes under tension. 

Finally, the effective linear spring model [221134] is a very useful simplification of membrane elastic 
models when dealing with local thickness deformations and hydrophobic mismatch. Its applicability has 
been thoroughly tested on systems where gramicidin is used to probe the influence of various molecules 
on membrane properties (see, e.g., Ref. |10j). As other quadratic elastic models, our model yields an 
effective spring model. However, since the expression of the spring constant depends on the details of the 
model, careful consideration is required when one is interested in the behavior of a particular material 
constant. 



Methods 

1 Derivation of the effective Hamiltonian 

1.1 General expression of the bilayer effective Hamiltonian 

Let us consider a patch of bilayer membrane with a fixed projected area A p , at fixed chemical potential 
fi. The rest of the membrane (e.g., of the vesicle) plays the part of the reservoir that sets the chemical 
potential fi. The effective Hamiltonian per unit projected area in each monolayer is = /^/S^, where 
is given by Eq. [TJ while the projected area S ± per molecule reads = E ± [l — (V/i ± ) 2 /2] to second 
order. Hence, Eq. Q] yields, to second order in the deformation and in the relative stretching of the 
monolayers, 

/± .^ + f So (5i^) 2 4^4(^) vV 

4 2^0 2^0 ^0 ^ ^0 

We assume that the hydrophobic chains of the lipids are incompressible. Let us introduce the excess 
hydrophobic thickness u + (resp. u~) of the upper (resp. lower) monolayer, defined as its hydrophobic 
thickness along the normal to its hydrophilic-hydrophobic interface minus the equilibrium monolayer 
hydrophobic thickness do/2 (see Fig. [T]). In the spirit of Refs. [T2TU4j . we use the incompressibility 
condition 

v = S± («± + |) , (21) 

where v is the constant hydrophobic volume per lipid. (In this incompressibility condition, a correction 
arising from membrane curvature is neglected. Using the complete incompressibility condition instead 
of this one yields the same effective Hamiltonian Eq. [2] but with different expressions of c$ and k as a 
function of the constants involved in Eq. [T] These expressions depend on /i, and consequently on the 
applied tension, but this dependence is negligible for realistic tension values. As the rest of our discussion 
is not affected by this, we keep the approximate incompressibility condition for the sake of simplicity. 
Note that the exact incompressibility condition was implemented recently in Ref. |23j.) 

In all the following, we will work to second order in the small dimensionless variables u /do, \Vu \, 
doV 2 u ± , |V/i ± | and do^ 2 ^. In this framework, using the relations (S — S*) /S ± = 2u ± /d and 
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m ± V 2 w ± = V (u ± Vm ± ) - (Vm ± ) 2 , Eg .[2U1 becomes 

2ti± (V/i±) 2 \ K 



r 



S 



do 



o 



I + ^( u ±) 2 ± ^V 2 ^ ± ^(c - c E ) u± V 2 ^ 



^(v 2 /i ± ) 2 + |dct(a J a J / l ± ) 



— V • (m ± Vu ± ) - V 2 m ± 
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+ «(V 2 U ± ) 2 



(22) 



In this expression, we have introduced the constitutive constants of a monolayer: /q Eo = K a /2 is 
compressibility modulus of the monolayer, /2/(2Eo) = «o/2 is its bending rigidity, /k/Eo = k/2 is 
its Gaussian bending rigidity, /1//2 = Co is its spontaneous (total) curvature, and /1//2 = c is the 
modification of the spontaneous (total) curvature due to area variations. More precisely, c' = dc s /dE 
where c s (E) = co + c (E — Eo) is the lipid area-dependent (total) spontaneous curvature of the monolayer. 
In addition, recall that do denotes the equilibrium hydrophobic thickness of the bilayer membrane. Finally, 
we have introduced the constants 



a E + (3 



dl 



k" = 4 



CSo 
dl 



(23) 
(24) 



These two constants have the dimension of a surface tension, like K a . 

In our description, the state of monolayer ± is determined by the two variables and u ± . Hence, 
the state of the bilayer membrane is a priori determined by four variables. However, given that there 
must be no space between the two monolayers, the distance along z between the hydrophilic-hydrophobic 
interfaces of the two monolayers must be equal to the sum of their projected thicknesses. Hence, to 
second order, we have the following geometrical constraint: 



h+ - h~ = 



do 
2 



1 - 



(V/i H 



do 
2 



1 - 



(25) 



This leaves us with only three independent variables to describe the state of the membrane. Let us choose 
the average shape h of the bilayer, the sum u of the excess hydrophobic thicknesses of the two monolayers, 
and the difference <5 between them: 



h = 



u = u + + u 



(26) 

(27) 
(28) 



Thus, we can rewrite the effective Hamiltonian / = /+ + /~ per unit projected area of the membrane 
in terms of the new variables h, u and 5. It reads, to second order in the small dimensionless variables 
u/do, 8 /do, |Vm|, |V5|, |V/i|, doV 2 u, doV 2 <5, and doV 2 h, and discarding derivatives of order higher than 
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two: 



f = o 



1 



d~ + 2 



8 



(V^) 2 + i(V 2 w) 2 



+ 2HJ± v 2 u + (c - c £ ) (« V 2 u + 2 5 V 2 ^) 



det(didjh) + - det(9 i 9 7 -u) 



£[(Vt 



(V^) 2 



fr" d 2 



[(V 



2 u) 2 



(V 2 <5) 2 1 



^ V 2 u + ^ [ V • («Vu) + V • (<SW) ] , 
do a 



(29) 



where we have introduced a 
Methods, Sec. [2). 



-2/i/So, which plays the part of an externally applied tension (see 



1.2 Eliminating <5 

In the present study, we are not interested in the variable 8. In a coarse-graining procedure, this degree of 
freedom can be eliminated by integrating over it. In our Gaussian theory, it simply amounts to minimizing 
/ with respect to 8. This variable is coupled to the membrane curvature V 2 /i, but not to u. In the case 
of a constant curvature, the constant value 



do kq . 2 



(30) 



is a simple solution to the Euler-Lagrange equations in 8, for which the term involving 8 in / reads 



fs = -\^(co-c'^o) 2 (V 2 /,) 2 



(31) 



As the variable 8 varies spontaneously on length scales much shorter than the variable h, we can consider 
in a first approximation that S will simply follow V h, in which case this constant solution is the valid 
one. Thus, after this partial minimization, this term provides a correction to no- 
Wc finally obtain 



f = <? 



u 
do 



(Vhf 



Ka 2 

2d% 



y + 8"2^ (c °- c °^ 



(V«) 2 



k , Kd\ 



(v 2 u y + n 



1 



det(didjh) + — det(didju) 



Ko Co 




V 2 u + 


2 


d _ 





— (co-c E ) + ^ 



V ■ (uVu) ■ 



(32) 



where the usual Helfrich bending rigidity k, associated with the average shape, is related to kq through 



/ (rr \! 

k = no - —(c - c L ) 

a. 



(33) 



In the case where the average shape of the membrane is flat, i.e., h = 0, dropping constant terms, we 
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obtain the expression of / in Eq. [5] with 

K = -^(co-c'^o) + k' a + ^, (34) 
K = ^- + Kdl, (35) 

Kq Cq 2/3 

Al = — -^' (36) 

Thus, in general, in Eq. [21 the constants K' a , K" include contributions in k' a and A;", which arise from 
a, f3 and ( (see Eqs. [23] [24]). Therefore, the terms in gradient and Laplacian of £ introduced in Eq. [T] 
cannot be neglected a priori, as they contribute to the terms in (Vit) 2 and (V 2 «) 2 that are traditionally 
accounted for in models describing membrane thickness deformations [7J [T^HHl HH1 ■ Due to these 
contributions, the values of the constants K' a and K" are not fully predicted by the constants involved 
in the Helfrich model. This stands in contrast with the models developed previously [71 [T2T[T^l[T51[2Tj . In 
addition, the terms arising from a, (3 and ( modify the relations between the various coefficients: in the 
previous models that accounted for boundary terms, assuming a = (3 = £ = 0, and disregarding tension, 
one had K' a = — 2 [14], which is no longer true here. This will affect the equilibrium thickness profile 
of a membrane containing a mismatched protein. 



1.3 Link with the Helfrich Hamiltonian 

Since the variables h and u are decoupled in the Hamiltonian density / given by Eq. 1321 the terms 
depending on h can be isolated, yielding 

+ \ (V 2 ^) 2 + Rdetididjh) , (38) 

which corresponds to the Helfrich Hamiltonian [5] for a membrane composed of two identical monolayers. 
In particular, the term in a has the standard form of a Helfrich tension term, conjugate to the actual area 
A of the membrane, since the element of area is dA = dxdy^Jl + (V/i) 2 = dxdy [l + (V/i) 2 /2] to second 
order. Hence, a can be viewed as an effective applied tension. This interpretation of a is explained in 
more detail in Sec. [2] of our Methods part. 

Hence, our model gives back the Helfrich Hamiltonian if the state of the membrane is described only 
by its average shape h, i.e., if the variable u is integrated out. 



fh 



(vhy 



1.4 Stability criterion 

Let us focus on a membrane with flat average shape h, described by Eq. [5] Depending on the values 
of the constants K a , K' a and K", a homogeneous thickness u — can be less or more energetically 
favorable than an undulated shape. The physical situation we wish to describe is the one where the 
equilibrium state has a homogeneous thickness. To determine which sets of constants comply with this, 
let us calculate the effective Hamiltonian per unit projected area /def of a membrane with harmonic 
undulations characterized by the wave vector q. Neglecting boundary terms (by taking appropriate 
boundary conditions or by assuming that the undulations decay on some large length scale), we obtain 
/dof oc K a I ' d\ + K' a q 2 + K"q A , where the omitted prefactor is positive. The flat shape is favored if /d c f > 
for all q, and otherwise there exist some values of q for which it is unstable. Thus, the conditions for the 
stability of the flat shape are K a > 0, K" > and K' a > -2^/ K a K[[ / d . 
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2 Membrane submitted to an external tension 

In Sec. [T] of our Methods part, we have derived the effective Hamiltonian of a bilayer membrane in the 
(//, A p ) ensemble. This is the most convenient thermodynamic ensemble to work in. However, in order 
to describe experiments where a vesicle is submitted to an external tension, one should work in the 
(iV, t) ensemble, where N is the number of lipids in the vesicle and r is the externally applied tension. 
This is especially interesting in order to analyze the results of Ref. [TB]. The ensemble change can be 
performed using a Legendre transformation: in the (iV, r) ensemble, the adapted effective Hamiltonian 
is G(N, t) = F(fx, A p ) + fiN — tA p , where F(/j,, A p ) = J A dxdy f, with / expressed in Eq. [32] and 



OjJL 



dF 



dA r , 



(39) 



Let us restrict ourselves to the case of a homogeneous and flat membrane, i.e., to a membrane with 
constant h and u. Then, using Eq. [S^to eliminate the variables /i and A p from the expression of G, we 
obtain, to second order: 



G(N,r) = N^- 
do 



U t T r „ S U 



(40) 



Minimizing G with respect to u yields the equilibrium excess thickness u eq of the membrane at a given 



imposed tension r. To first order, it reads 



T 

Xeq = ~—do, (41) 



Note that, since u/do is assumed to be a first-order quantity, r/K a must be first-order too for our 
description to be valid for u = u cq . This property has been used to simplify the result in Eq. 1411 In 
practice, r <C K a is well verified, given that r cannot exceed a few mN/m without the vesicle bursting, 
while K a is of order 100 mN/m. Since do is the equilibrium hydrophobic thickness of this piece of 
homogeneous and flat membrane submitted to a vanishing external tension, it is consistent that u eq 
vanishes when r does, as u is the excess thickness with respect to do- Eq.@J] shows that the thickness of a 
membrane with fixed number of lipids decreases when the external tension increases, and is in agreement 
with Ref. PI]. 

We are now going to show that the constant a in the (n, A p ) ensemble (see, e.g., Eq. 13"!?]) plays the 
part of an externally applied tension. For this, let us calculate the equilibrium thickness of a membrane 
patch with projected area A p at a chemical potential \x, when it is homogeneous and flat. This amounts 
to minimizing / with respect to u. For a homogeneous and flat membrane, Eq. 1321 becomes 



u \ K„ u 

'='l 1+ 4 + iV (42) 







Minimizing / with respect to u then gives 



a 



, d . (43) 



Comparing Eq. [43] to Eq. [41] shows that a plays the part of the externally applied tension r. Hence, a 
can be considered as an effective applied tension. 

3 Membrane containing a cylindrical mismatched protein 

In this Section, we write down explicitly the equilibrium shape and the deformation energy of a membrane 
which contains a single cylindrical transmembrane protein with a hydrophobic mismatch (see Fig. [TJ3) . 



20 



This protein can correspond to a gramicidin channel in the dimer state. We focus on a membrane with 
a fiat average shape, described by the effective Hamiltonian per unit projected area in Eq. [2J We denote 
the radius of the protein by tq, and its hydrophobic thickness by I. We take the center of the cylindrical 
protein as the origin of the frame, which yields cylindrical symmetry. 

In order to treat the case where the membrane is submitted to a tension a, we rewrite Eq. [2] in terms 
of the variable u = u — u cq = u + ado/K a , which represents the excess hydrophobic thickness of the 
bilayer relative to its equilibrium value at an applied tension a (see Eq. 143 j). Discarding constant terms 
and using the relation a <C K a , which yields ad^A^i 'K a <C A\, it yields 

/ = l|" 2 + lf (v-) 2 + ^(v 2 ^) 2 

+ Ai V 2 ii + A 2 V • (uVu) + ^ det^dj-u) . (44) 



3.1 Equilibrium thickness profile 



Let us first review (see, e.g., Ref. |22j ) the equilibrium thickness profile u of the membrane containing 
the mismatched protein. This equilibrium shape is solution to the Euler-Lagrange equation associated 
with the effective Hamiltonian in Eq. |4"41 



V 4 m 



K" 



K a _ 

u = . 



K"d 2 



(45) 



Using the cylindrical symmetry of the problem and choosing solutions that vanish at infinity, we 
obtain, if the stability condition Eq. |8] is verified, the following solution to the Euler-Lagrange equation 
Eq.EU 

2(r) = A + K (k+r) + A_K {k^r) , (46) 
where K„ is the n th -order modified Bcssel function of the second kind, and 




K 

K" 



K a d l 



1/2' 



1/2 



(47) 



which are either both real or complex conjugate. 

The integration constants A± are determined by the boundary conditions at r = tq. The first bound- 
ary condition corresponds to strong hydrophobic coupling: on the inclusion boundary, the hydrophobic 
thickness of the membrane is equal to that of the inclusion, which is denoted by £ (see Fig.[T)3). It yields 
u(tq) = uq = £ — do (to first order, as explained in our Section entitled "Deformation profiles close to a 
mismatched protein"), or equivalently u(j"o) = uq = £ — do (1 — a / K a ). As far as the second boundary 
condition at r = ro is concerned, we will treat explicitly two different cases, which correspond respectively 
to a fixed slope and to a free slope in tq, as explained in the main text of the article. 



Fixed slope. In the case where the boundary conditions in r = are 

KJ , (48) 



u(r ) = u = £ - d ( 1 



u (r J = s 

which corresponds to a strong hydrophobic coupling and a fixed slope s at r = tq, we obtain: 

K^s + k T Kfu 



T • 



(49) 
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whore 

K± = K„ (fc±r„) . (50) 

Note that A + and A_ are either both real or complex conjugate (like k±), which ensures that the solution 
Eq. [46] is real. 

Free slope. An alternative choice of boundary conditions in r = tq is 

(51) 



to first order again. The first of these conditions corresponds to a strong hydrophobic coupling, as before. 
The second one arises from minimizing the total free energy of the system without further constraints. It 
corresponds to the case where the slope at r — ?*o is free to adjust itself to yield the smallest deformation 
energy. With these "free-slope" boundary conditions, we obtain: 

R MCfjlo - 4r K T [Ai + (A 2 + K>'k\) u ] 
± 4r i^'(fc2 _fc2)K+Ko -S(fc+KoK+-fc_K+Kr) ' 1 ' 

which arc, again, cither both real or complex conjugate. 

Let us now assume that ft — £ = 0, as in the main text of this article. In order to understand the 
impact of k' a (i.e., of a) on A± in the free-slope case, let us express A± as a function of k' a , ro, do and 
of the bulk constants K a , K' a and K'^, whose values can be extracted from the fluctuation spectra in 
simulations. Using Eq. E7J the relation A\ = 2K"cq, which can be derived from Eqs. [7] and and the 
relation A 2 = (k' a — K' a )/2, which stems from Eqs. 2] and [71 we obtain: 

A ± Rk T Kfu - 2r KT {AK'^cq + [k' a ± K'j (kl - k\)] u Q ) 
± 4r K» (k 2 + - k 2 _ ) K+Kq - R (k + KoK+ - k-K+K^) 

For fixed values of ro, do, K a , K' a and K'^, the constants A± can be viewed simply as functions of k' a and 
Co: let us denote them by A±(k' a , cq). The following relation holds for all k' a and Co: 

A±{k' a ,co) = A±(0,co) , (54) 

with 

k' 

co = c + T^77«o ■ (55) 

Hence, in the framework of a model that assumes k' a = 0, the effect of a nonvanishing k' a on the equilibrium 
membrane thickness profile would be that Co is replaced by a renormalized spontaneous curvature Co, 
which depends linearly on uq. At vanishing applied tension (in which case, uq = uq), and neglecting the 
difference between n = K"/4 and k, we obtain Eq. fl6l 



3.2 Deformation energy 

Let us now calculate the deformation energy F of the membrane due to the presence of the mismatched 
protein. For the equilibrium shape of the membrane, which is solution to the Euler-Lagrange equation 
Eq. I45[ we are left only with boundary terms at the inclusion edge in r = ro (no other boundary terms 
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contribute, since the deformation u caused by the presence of the mismatched channel vanishes sufficiently 
far away from it). We can write 



F = 



dxdy f = 2n / rdr f 

I r 



ttI K"r 



dr 



V 2 u) -u'V 2 u 



A\ r u + A 2 r u u + — u' 2 
8 



K 

K" 



r=r 



(56) 



where 



du/dr. We have used the expression of the Gaussian curvature for small deformations in a 



system with cylindrical symmetry: dct(didjit) = u'u" jr = (2r) 1 d(u' 2 )/dr. To express the deformation 
energy F explicitly, one has to use the boundary conditions in r = ro- 

Fixed slope. For the boundary conditions in Eq. 051 corresponding to a fixed slope in ro, using 
Eqs. 1461 1471 and 1491 we can rewrite the deformation energy of the membrane in Eg. 1561 as 



-2tt 



Ax r s + A 2 r u s + - s 



/c+Kq K+ - fc-K^Ki 
2/fc + fc_ (k+K+K- - fc_Kp K+) M 



k + k_ (k\ - k 2 _) K+K^ 

K+K- (k 2 + - k 2 _) s 2 



(57) 



This expression shows that F is a second-order polynomial in uq and s. 

Spring constant for s = 0. In the particular case where the fixed slope s vanishes, Eq. [57] becomes 



F = H a u 2 



o « ' 



where the effective spring constant reads 

nr K" k+k^ (k\ - k 2 _) K+Kf 



Hn = 



(58) 



(59) 



Dependence on applied tension. Since uq = I — do(l — <j/K a ), Eq. [571 shows that F is a second- 
order polynomial in the applied tension a. (In our model, K' a features a contribution coming from a, see 
Eq. 2J However, as mentioned in the main text, the dependence of K' a on a is negligible in practice, and 
we thus disregard it: in this framework, C\ and C 2 do not depend on a.) We can write 



F 
kWT 



= C + C x o + C 2 a 2 



(60) 



with 



C x = 



2wd r K" k + k- 



k B TK a (fc_K+Kf 



/c + K -K+) 



{k+K+K^ - fc_Ko Kf) s 



{k\ - k 2 _) KrK+ (d - I) 



2nd r Q 
k B TK a 



sAo 



G 2 = - 



d 2 . Hp 
K 2 knT ' 



(61) 
(62) 
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where Hq is the effective spring constant expressed in Eq. 1591 Note that k and A\ do not appear in the 
coefficients C\ and C 2 , and that A2 and s are oniy present in C\. 

Free slope. For the boundary conditions in Eq. 1511 corresponding to a free slope in ro, using 
Eqs. l46l l47l and [52l we can rewrite the deformation energy of the membrane (see Eq. [56]) as 



F 



nr 



K (k + K^K+ - fc-K+Kr) - 4r K» (k\ - k 2 _) K+K Q - 
4r (fc+K -K+ (A 2 + K'^k 2 _f - fc^K+Kf (A 2 + K"K 



K'^R (k\ - k 2 _) k+k-KfK^ 



8A iro 



A 2 (fc+KoK+ - fc-K+Kf) u + Ul\tq {k+^K+ - fc-K+K^) 



(63) 



This expression shows that F is a second-order polynomial in uq. 
Spring constant. Eq. [63] can be expressed as 

F = H f (u Q -uf n ) 2 + F mh 

where the effective spring constant reads 

7rr 



H 



I 



R(k+K^K+ - fc-K+K^) 



kl) K+K 



4r K'> (k 2 

4r (fc+K-K+ (A 2 + K'^kl) 2 - fc-K+K^ (A 2 + K'^k 
f K'^K {k\ - k 2 _) fc+fc_K+K7 



(64) 



(65) 



while u™ ln denotes the value of uq that minimizes F, and F min is the minimum of F, obtained for uq = 
u™ 111 . Note that both it™ 111 and F mm are nonzero if A\ ^ (see Eq. lSB"]) , due to the spontaneous curvature 
of each monolayer. The effect of monolayer spontaneous curvature was disregarded in Ref. [T8ll22| . which 
explains why Eq. [64] differs from the standard expression F = HfU^ |22] , 

Dependence on applied tension. Since uq = (. — do(l — <r/K a ), Eq . 1551 shows that F is a second- 
order polynomial in the applied tension a (neglecting the cr-dependence of K' a as explained in the main 
text). Thus, we can write 

- -J-= = Co + C x u + C 2 a 2 , (66) 
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with 



Ci 



-2irr a d 

k B TK a [k (fc + K"K+ - fc_K+K") - 4r K» (fc 2 - k 2 _) K+K ] 



4A iro 



A 2 (fc + K -K+ - fc_K+Kr) + 4r (fc+K -K+(A 2 +^'fc 



fc_K+Kf (A 2 + + K%k (kl - k 2 _) k+k-K+K^ {I - do) 



Co 



2 a n f 



(67) 
(68) 



Kl k B T ' 

where Hf is the effective spring constant expressed in Eq. 1651 
4 Estimating c' 

Let us start from the free energy per molecule in monolayer + expressed in Eq. [T] All the quantities 
involved in this expression are defined on the hydrophilic-hydrophobic interface S of the monolayer. 

Let us consider a surface S' parallel to S, and let us call 5 the algebraic distance from S' to S. To 
second order in the small dimcnsionless variables c\8 and c 2 <5, where c\ and c 2 denote the local principal 
curvatures of the monolayer (recall that H = (ci + c 2 )/2 and K = cic 2 ), geometry gives [19] : 



E' = £ (1 + 2 HS + KS 2 ) , 
H' = H + (K -2H 2 ) S , 
K ' = K . 

Hence, we can rewrite / + using variables defined on S 1 , to second order: 

/+ = - So) 2 + fiH' + (f{ - 2 f» E S) (E' - E ) H' 

+ (/ a + 2/^' E^ «5 2 - 2 /( E S + 2 h 5) H' 2 + (f K - h 5) K> 
+ a (VE') 2 + & V 2 E' + C (V 2 E') 2 - (i, , 



(69) 
(70) 
(71) 



(72) 



where we have neglected terms containing derivatives of order higher than two. 

If S' is the neutral surface of the monolayer [19] , by definition, the curvature and the area variations 
are decoupled, which entails f[ = 2 f' Q ' Eo £, where £ denotes the algebraic distance from the neutral 
surface to the hydrophilic-hydrophobic interface of the monolayer. Thus, given that /g' = K a /(2'So), 
f-2 = no Eq, and /{// 2 = cJ (see Methods, Sec. II. ip . we obtain 



c n E 



Kq 



(73) 



Acknowledgments 

We thank Grace Brannigan for sharing with us some data corresponding to the simulations described 
in Ref. [15]. We thank Mark Goulian for sharing with us a notebook containing some of the original 
calculations of Ref. [T5], and for email correspondence. We also acknowledge critical reading of our 
manuscript by Florent Bories. 



25 



References 

1. Mouritsen OG (2005) Life - as a matter of fat: the emerging science of lipidomics. Springer, Berlin. 

2. Sackmann E (1984) Physical basis for trigger processes and membrane structures. In: Chapman 
D, editor, Biological Membranes, London: Academic Press, volume 5. pp. 105-143. 

3. Hclfrich W (1973) Elastic properties of lipid bilaycrs - Theory and possible experiments. Zeitschrift 
fur Naturforschung C - Journal of Biosciences 28: 693-703. 

4. Owicki JC, Springgate MW, McConnell H (1978) Theoretical study of protein-lipid interactions in 
bilayer membranes. Proc Natl Acad Sci USA 75: 1616-1619. 

5. Owicki JC, McConnell H (1979) Theory of protein-lipid and protein-protein interactions in bilayer 
membranes. Proc Natl Acad Sci USA 76: 4750-4754. 

6. Mouritsen OG, Bloom M (1984) Mattress model of lipid-protein interactions in membranes. Bio- 
phys J 46: 141-153. 

7. Huang HW (1986) Deformation free energy of bilayer membrane and its effect on gramicidin channel 
lifetime. Biophys J 50: 1061-1070. 

8. Kelkar DA, Chattopadhyay A (2007) The gramicidin ion channel: A model membrane protein. 
Biochim Biophys Acta - Biomembr 1768: 2011-2025. 

9. O'Connell AM, Koeppe II RE, Andersen OS (1990) Kinetics of gramicidin channel formation in 
lipid bilayers: transmembrane monomer association. Science 250: 1256-1259. 

10. Lundbaek JA, Koeppe II RE, Andersen OS (2010) Amphiphilc regulation of ion channel function 
by changes in the bilayer spring constant. Proc Natl Acad Sci USA 107: 15427-15430. 

11. Hclfrich P, Jakobsson E (1990) Calculation of deformation energies and conformations in lipid 
membranes containing gramicidin channels. Biophys J 57: 1075-1084. 

12. Dan N, Pincus P, Safran SA (1993) Membrane-induced interactions between inclusions. Langmuir 
9: 2768-2771. 

13. Aranda-Espinoza H, Berman A, Dan N, Pincus P, Safran S (1996) Interaction between inclusions 
embedded in membranes. Biophys J 71: 648. 

14. Brannigan G, Brown FLH (2006) A consistent model for thermal fluctuations and protein-induced 
deformations in lipid bilayers. Biophys J 90: 1501-1520. 

15. Brannigan G, Brown FLH (2007) Contributions of gaussian curvature and nonconstant lipid volume 
to protein deformation of lipid bilayers. Biophys J 92: 864-876. 

16. West B, Brown FLH, Schmid F (2009) Membrane-protein interactions in a generic coarse-grained 
model for lipid bilaycrs. Biophys J 96: 101-115. 

17. Elliott JR, Needham D, Dilger JP, Haydon DA (1983) The effects of bilayer thickness and tension 
on gramicidin single-channel lifetime. Biochim Biophys Acta - Biomembr 735: 95-103. 

18. Goulian M, Mesquita ON, Fygenson DK, Nielsen C, Andersen OS, et al. (1998) Gramicidin channel 
kinetics under tension. Biophys J 74: 328-337. 

19. Safran SA (1994) Statistical thermodynamics of surfaces, interfaces and membranes. Addison- 
Wesley. 



2G 



20. Bitbol AF, Peliti L, Fournier JB (2011) Membrane stress tensor in the presence of lipid density 
and composition inhomogencities. Eur Phys J E 34: 53. 

21. Fournier JB (1999) Microscopic membrane elasticity and interactions among membrane inclusions: 
interplay between the shape, dilation, tilt and tilt-difference modes. Eur Phys J B 11: 261-272. 

22. Nielsen C, Goulian M, Andersen OS (1998) Energetics of inclusion-induced bilayer deformations. 
Biophys J 74: 1966-1983. 

23. Watson MC, Pencv ES, Welch PM, Brown FLH (2011) Thermal fluctuations in shape, thickness, 
and molecular orientation in lipid bilayers. J Chem Phys 135: 244701. 

24. Israclachvili JN (1992) Intcrmolccular and surface forces, Second edition. Academic Press. 

25. Sharp KA, Nicholls A, Fine RF, Honig B (1991) Reconciling the magnitude of the microscopic and 
macroscopic hydrophobic effects. Science 252: 106-109. 

26. May S, Ben-Shaul A (1999) Molecular theory of lipid-protein interaction and the L a -H| transition. 
Biophys J 76: 751-767. 

27. Hladky SB, Gruen DWR (1982) Thickness fluctuations in black lipid membranes. Biophys J 38: 
251-258. 

28. Rawicz W, Olbrich KC, Mcintosh T, Needham D, Evans E (2000) Effect of chain length and 
unsaturation on elasticity of lipid bilayers. Biophys J 79: 328-339. 

29. May ER, Narang A, Kopelevich DI (2007) Role of molecular tilt in thermal fluctuations of lipid 
membranes. Phys Rev E 76: 021913. 

30. May ER, Narang A, Kopelevich DI (2007) Molecular modeling of key elastic properties for inho- 
mogeneous lipid bilayers. Mol Simulat 33: 787-797. 

31. Watson MC, Brandt EG, Welch PM, Brown FLH (2012) Determining biomcmbrane bending rigidi- 
ties from simulations of modest size. Phys Rev Lett 109: 028102. 

32. Lindahl E, Edholm O (2000) Mesoscopic undulations and thickness fluctuations in lipid bilayers 
from molecular dynamics simulations. Biophys J 79: 426-433. 

33. Marrink SJ, Mark AE (2001) Effect of undulations on surface tension in simulated bilayers. J Phys 
Chem B 105: 6122-6127. 

34. Lundback JA, Andersen OS (1999) Spring constants for channel-induced lipid bilayer deformations 
estimates using gramicidin channels. Biophys J 76: 889-895. 

35. Kim T, Lee KI, Morris P, Pastor RW, Andersen OS, et al. (2012) Influence of hydrophobic mismatch 
on structures and dynamics of gramicidin A and lipid bilayers. Biophys J 102: 1551-1560. 

36. Kolb HA, Bamberg E (1977) Influence of membrane thickness and ion concentration on the prop- 
erties of the gramicidin A channel: autocorrelation, spectral power density, relaxation and single- 
channel studies. Biochim Biophys Acta - Biomembr 464: 127-141. 

37. Hwang TC, Koeppe II RE, Andersen OS (2003) Genistein can modulate channel function by a 
phosphorylation-indcpcndcnt mechanism: importance of hydrophobic mismatch and bilayer me- 
chanics. Biochemistry 42: 13646-13658. 



27 



38. Chung H, Caffrcy M (1994) The curvature elastic-energy function of the lipid- water cubic 
mesophase. Nature 368: 224-226. 

39. Chung H, Caffrey M (1994) The neutral area surface of the cubic mesophase: location and prop- 
erties. Biophys J 66: 377-381. 

40. Tcmpler RH (1995) On the area neutral surface of inverse bicontinuous cubic phases of lyotropic 
liquid crystals. Langmuir 11: 334-340. 

41. Templer RH, Turner DC, Harper P, Scddon JM (1995) Corrections to some models of the curvature 
elastic energy of inverse bicontinuous cubic phases. J Phys II France 5: 1053-1065. 

42. Vacklin H, Khoo BJ, Madan KH, Seddon JM, Tcmpler RH (2000) The bending elasticity of 1- 
monoolcin upon relief of packing stress. Langmuir 16: 4741-4748. 

43. Szule JA, Fuller NL, Rand RP (2002) The effects of acyl chain length and saturation of diacylglyc- 
erols and phosphatidylcholines on membrane monolayer curvature. Biophys J 83: 977-984. 

44. Leikin S, Kozlov MM, Fuller NL, Rand RP (1996) Measured effects of diacylglycerol on structural 
and elastic properties of phospholipid membranes. Biophys J 71: 2623-2632. 

Figure Legends 








h ■ ' 












h 




d 




I 














h . 












B 



Figure 1. Definitions. A) Cut of a bilayer membrane. The solid black lines mark the boundaries of 
the hydrophobic part of the membrane, and the exterior, which is shaded in blue, corresponds to the 
hydrophilic lipid heads and the water surrounding the membrane. The hydrophobic thickness, defined 
along the normal to the hydrophobic-hydrophilic interface, of the upper (resp. lower) monolayer, shaded 
in orange (resp. yellow), is u + + do/2 (resp. u~ + do/2). The height of monolayer ± along z is denoted 
by h ± . The average membrane shape, h = (h + + h )/2, is represented as a red dashed line. B) Cut of 
a bilayer membrane (with hydrophobic part shaded in yellow) containing a protein with a hydrophobic 
mismatch (orange square). The equilibrium hydrophobic thickness of the bilayer is do, while the 
hydrophobic thickness of the protein is £. The average shape of the membrane is flat, and the thickness 
deformations of the two monolayers are identical (u + = u~ = u/2). Hence, the average shape h is 
constant, and confounded with the midlayer of the membrane. Although u ± is defined along the normal 
to the monolayer hydrophilic-hydrophobic interface, the boundary condition at the inclusion edge, i.e., 
in r = ?*o, simply reads u(tq) = uo = £ — do to first order (see main text, Section entitled "Deformation 
profiles close to a mismatched protein" ) . 
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Figure 2. Thickness gradient. Cut of a bilayer membrane with a symmetric thickness gradient. The 
dashed blue lines correspond to the hydrocarbon- water interfaces. 
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Figure 3. Thickness deformation due to a mismatched inclusion. Membrane thickness profile 
from Ref. |16| in the vicinity of a mismatched inclusion with hydrophobic thickness £ ~ 2.4 nm and 
radius tq = 9 A, with center in r = 0, as a function of the radial coordinate r. The equilibrium 
membrane hydrophobic thickness is do — 3.6 nm. The unit of length on the graph is 6 A, as in Ref. [IB] 
Dots: numerical data (the error bars on the data, not reproduced here, are about 1 A wide |16jV Line: 
best fit. Exactly as in the original reference, the numerical data is fitted to Eqs. I46ll53l with k' a = 0, 
taking uq and the (renormalized) spontaneous curvature cq as fitting parameters, the other constants 
being known from the fluctuation spectra. 
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Figure 4. Renormalized spontaneous curvature cq as a function of the hydrophobic 
mismatch uq. Data from Ref. [16] . which presents fits of simulation results for inclusions with three 
different hydrophobic thicknesses. Line: linear fit, with slope (0.56 ± 0.02) nm~ 2 . Note that our Co 
corresponds to twice that in Table 2 of Ref. [16] , as we work with total curvatures instead of average 
curvatures. The error bars on c are those listed in that table, and uq corresponds to 2 in that table. 
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Figure 5. Renormalized spontaneous curvature cq as a function of the inclusion radius r$ 
and the hydrophobic mismatch uq. A) Co versus rg. The values of Co were obtained by fitting each 
thickness deformation profile of Ref. [15] . Circles (blue): positive mismatch, «o = 0.45 nm. Squares 
(red): negative mismatch, uo = — 1.1 nm. Solid lines: average values; dotted lines: standard deviation 
over the seven data points (corresponding to the different ro), for each value of uo- B) Average value of 
Co (see A) as a function of the hydrophobic mismatch uq. The equation of the line joining the two data 
points has a slope 0.26 nm - 2 . 
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Figure 6. Renormalized spontaneous curvature cq versus the relative volume variation 

v(ro)/vo on the inclusion edge. The values of Co are extracted from fitting the data of Ref. [15], and 
the values of v(tq)/vq are directly taken from Ref. [15] . 
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Figure 7. Formation rate / of gramicidin channels versus the applied tension a, analyzed 
with a quadratic model. Diamonds: experimental data retrieved from Fig. 6b of Ref. [18] . after 
subtraction of the baselines Cq . Dotted black line: best quadratic fit, with C\ = 0.74 (mN/m)" and 
Ci = —0.09 (mN/m)~ ; y 2 = Xmin- Dashed green line: results obtained from the elastic model of 
Ref. (22], with the constants given in [183; X 2 /Xmin = 5-72. Dashed-dotted blue line: idem with more 
recent values of the constants; X 2 /Xmin = 6.68. Solid red line: results obtained by taking s = and the 
recent values of the constants in the model of Refs. 0(23 ; X 2 /Xmin = The values of C\ and Ci 

corresponding to the curves on this graph are listed in Table [2] 
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Figure 8. Comparison between the experimental values of C\ and C-2 and those obtained 
from different models. Colorscale: goodness-of-fit function \ 2 (see Eg. [T9|) for the data of Ref. |18j . 
as a function of the fitting parameters C\ and Ci. White diamond: values of C\ and C2 that give the 
best fit. Black triangle: results obtained from the elastic model of Ref. [22], with the constants given 
in [T5]. Lines: trajectories obtained from our model in the (Ci,C2) plane when varying k' a . Red: free 
slope; green: s = 0, black: s = 0.3. These three curves start by a white dot at k' a = 0, and k' a increases 
rightwards along these curves. The rightmost white dot (k' a = 0, s = 0) roughly corresponds to the best 
agreement we can obtain between our model and the experiment fitted to the quadratic model (red 
curve on Fig. [7]). The black diamond corresponds to the best agreement we can obtain between our 
model and the experiment fitted to the linear model at low tension (see Fig. [9|) . 



-0.2 



0.0 0.5 1.0 1.5 2.0 

o[mN/m] 

Figure 9. Formation rate / of gramicidin channels as a function of the applied tension a, 
analyzed with a linear model, for a < 2mN/m. Diamonds: experimental data retrieved from 
Fig. 6b of Rcf. [IS], after subtraction of the baselines Cq (which arc different from those of Fig. [7] since 
the fitting model is here linear instead of quadratic). Line: best linear fit, yielding 
Ci = (0.62 ± 0.05) (mN/m)" 1 ; correlation coefficient: r = 0.894. 
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Tables 







Set 1 


Set 2 


Set 3 


Free s 


fl> if k' a = 


41 


46 


33 


Free s 


A£ if ff/ = i/ oxp = 115 


25 


24 


26 


s = 


i/0 if = 


130 


133 


91 


s = 


k' a if Ho = H cxp = 115 


< 


< 


7.5 



Table 1. Spring constant H and constant k' a of monoolein. The results are given both for the 
free-slope boundary condition (using Eq. I65|) and for the zero-slope boundary condition s = (using 
Eg. I5"9"]) . All values of H and k' a are given in mN/m. Negative values of k' a are not detailed since they 
would yield an instability for the monolayer Hamiltonian Eq. [221 in the present framework where /c" = 0. 
The different columns correspond to three different data sets for the parameters of the membrane. Set 1 
corresponds to the data from [38] at 25°C: k = 3.6 x 10 -20 J, c = -0.13511m- 1 , R = 8.8 x 10" 22 J. Set 
2 takes into account the corrections on Co and R in [41]: Co = — 0.196 nm , R = —3.6 x 10~ 21 J. Set 3 
corresponds to the data from [32]: k = 1.2 X 10~ 20 J, c = -0.503 nm" 1 , and R = —1.2 x 10" 21 J 
deduced from R/k = —0.1 [41] . In all cases, we have taken tq = 1 nm |34j . do = 2.46 nm [35], 
C = -0.3A gD], H Q = 140mN/m [27][33|- 



Model 


Ref. [22] 


Ref. [21 


Ref. [21 


Ours, 
with 

K = o 


Ours, 
with 
^ = 


Ours, 
with 

K = o 


Constants 


Ref. [TH] 


Recent 


Recent 


Recent 


Recent 


Recent 


Slope s 


0.3 


0.3 





0.3 





Free 


Gi [10- 3 (mN/m)" x ] 


354 


282 


480 


292 


502 


339 


-C a [10- 3 (mN/m)" z ] 


21.4 


6.11 


6.11 


6.40 


6.40 


3.34 




5.72 


7.15 


1.75 


6.68 


1.75 


4.31 



Table 2. Values of C\ and C*2 obtained from the model of Ref. [22j and from our model 
with k' a = 0. The results are presented both for the fixed-slope boundary condition (see Eqs. ISTI 
and 162]) . with slopes and 0.3, and for the free-slope boundary condition (see Eqs. l67l and [68]). The 
corresponding values of % 2 are also given. Recall that the best quadratic fit to the data of Ref. [T5] 
yields d = 740 x lO" 3 (mN/m) _1 and C* 2 = -90.0 x 10~ 3 (mN/m)" 2 (see Fig. 0). 
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Slope s 





0.3 


Free 


k' a (mN/m) 





45 


30 


d [10- 3 (mN/m) L ) 


502 


490 


490 


-C 2 [10- 3 (mN/m)" 2 ] 


6.40 


9.17 


5.29 


x 2 


1.75 


1.69 


1.75 



Table 3. Values of k' a , C\ and C2 obtained from our model that yield the best agreement 
with the experimental results of Ref. |18j . analyzed with a quadratic fit (see Eq. 1181 and 
Fig. [7]). Results are presented for the fixed-slope boundary condition (see Eqs. I5T1 and I52"j) . with slopes 
and 0.3, and for the free-slope boundary condition (see Eqs. [67l and [68|) . 



Slope s 





0.3 


-0.17 


Free 


k' a (mN/m) 


23 


78 





60 


-C 2 [10- 3 (mN/m)^] 


7.90 


11.0 


6.39 


7.04 



Table 4. Values of k' a and C2 obtained from our model that yield the best agreement with 
the experimental results of Ref. [18j analyzed with the low-tension linear fit. More precisely, 
these values of k' a and Ci are associated with C\ = 0.62 (mN/m) - . Results are presented for the 
fixed-slope boundary condition (see Eqs. IBT1 and IB2"]) . with slopes 0, 0.3, —0.17, and for the free-slope 
boundary condition (see Eqs. l67l and 168)) . 



